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Abstract 

Purpose: This paper presents the preliminary results of a semi-automatic 
method for prostate segmentation of Magnetic Resonance Images (MRI) 
which aims to be incorporated in a navigation system for prostate brachyther- 
apy. 

Methods: The method is based on the registration of an anatomical 
atlas computed from a population of 18 MRI exams onto a patient im- 
age. An hybrid registration framework which couples an intensity-based 
registration with a robust point-matching algorithm is used for both atlas 
building and atlas registration. 

Results: The method has been validated on the same dataset that the 
one used to construct the atlas using the leave-one-out method. Results 
gives a mean error of 3.39 mm and a standard deviation of 1.95 mm with 
respect to expert segmentations. 

Conclusions: We think that this segmentation tool may be a very 
valuable help to the clinician for routine quantitative image exploitation. 

Key words: Atlas-based-Segmentation, prostate, MRI 

1 Introduction 

Adenocarcinoma of the prostate is the most common cancer expected to occur 
in men in 2008 1 . A total of 186,320 new cancer cases and 28,660 deaths from 
prostate cancer are projected to occur in the United States in 2008 pQ. 

Magnetic Resonance Imaging (MRI) allows the detection of the prostate 
and structures therein with higher accuracy than ultrasound images (US). Au- 
tomatic segmentations of the prostate in MRI have shown to be better corre- 
lated to manual expert segmentations than automatic segmentations in US. In 
order to make dose planning easier in prostate brachytherapy, we introduced 
a MRI/US registration method based on contours ([2] [3]); in the current ver- 
sion, segmentations are produced manually which is a practical limitation. The 



1 Excluding basal and squamous cell skin cancers and in situ carcinoma except urinary 
bladder 
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Figure 1: Overview of the method 
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objective is to automate this segmentation phase; in a first stage we focused 
on MM data. For practical use in a clinical situation, the algorithm must be 
accurate, relatively fast, and must be limited to as simple interaction with the 
expert as possible. 

Works on automatic or semi-automatic techniques for segmentation of the 
prostate in MRI are limited. Zwiggelaar and al. |I| have developed a technique 
based on polar-transform. Zhu and al. used an Active Shape Model to achieve 
a semi-automatic segmentation. Zhu and al. improve their previous work by 
using a combined 2D/3D Active Shape Model [5] [6]. [7] achieve an automatic 
segmentation using an atlas-based method. Other works deal with CT or US 
modalities : [5] use a coupled bladder and prostate segmentation in CT images 
to achieve an automatic prostate segmentation. [5] use an 2D active shape model 
to segment prostate in Ultrasound images. Our proposed segmentation method 
is based on the registration of an anatomical atlas to an individual image. The 
original aspect of this work consists to use a hybrid registration algorithm for 
atlas building and atlas matching. 



2 Material 

MRI acquisitions have been realized with a trans-rectal probe. Volumes have 
been acquired with a voxel size of 0.53x0. 53x3. 12 mm 3 and have been resamplcd 
to a voxel size of 1.0 x 1.0 x 1.0 mm 3 . The set of data used for model building 
is presently composed of 18 MRI exams. All acquisitions have been realized 
on patients where a prostate cancer was diagnosed, before a brachytherapy 
treatment. One expert has segmented all the MRI exams. These segmentations 
are used for model building and evaluation. 



3 Method Overview 

The proposed segmentation method is based on the registration of an anatom- 
ical atlas to an individual image. We both use an hybrid registration method 
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for atlas matching and atlas construction. A shape model of the prostate con- 
structed from the population used to build the atlas is used to regularize seg- 
mentations obtained after the atlas matching. Figure [T] illustrates the principle 
of the method. 

The hybrid registration which exploits image and geometric information is 
based on two alternated registrations. The framework registers a template image 
T(x) onto a reference image R(x) and a set of points M = {mo, m^} (called 
"model points" and belonging to the template image) to a set of points S = 
{s , sm} (called "scene points" and belonging to the reference image). The 
description of the general framework is done in section 14.11 Figure [2] illustrates 
the principle of the registration. 

During the anatomical atlas construction described in section 15.11 we use 
the hybrid registration framework to map all individuals of the database onto 
a common reference by exploiting an image similarity measure and a distance 
between expert segmentations. 

The same framework is also used during the atlas-based segmentation to 
incorporate a user interaction through a robust point to surface matching (see 
section UTTj) . The interaction is realized by asking the user to select one point on 
two specified regions in the image to be segmented. These regions are the base 
and the apex of the prostate and are defined by assigning a prior probability of 
match to each model point embedded in the atlas. This is described in section 

1531 

After the atlas-based segmentation a projection of the resulting segmentation 
on the shape space is done to regularize results (section [5J2]). A brief description 
of shape model construction is given in section RT^l 



4 An Hybrid Registration method 
4.1 A dual energy minimization 




Figure 2: Hybrid registration; left: principle of hybrid registration; right: 
correspondences between model points (circles) and a scene point (square) are 
modelled by assigning weights Wij to each couple (m,, Sj). 

Coupling intensity and geometric information can be very interesting to im- 
prove the quality of registration algorithms. Area of applications are for example 
the incorporation of segmentations extracted by some image processing tech- 
niques into an intensity-based registration ([10]). or the introduction of expert 
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constraints to an intensity-based registration ([11]). 

Here we propose a method that alternates the minimization of two criteria 
to both take into account an intensity-based distance and a geometric distance. 
The minimization of the intensity-based distance aims at matching the template 
image T(h(x)) on the reference image R{x) while the minimization of the ge- 
ometric distance aims at matching model points rrii belonging to the template 
image to scene points sj belonging to the reference images. The two criteria are 
incorporated into two dual energies functions which are minimized alternatively 
through two registrations algorithms: 

Eh. = E sim (T, R, h 1 ) (1) 
E 2 = ^^ U! . iJ ||/^ 1 (m l )-.s J || 2 +/3||/ l / -/i|| 2 (2) 

3 i 

Where h 1 and h are two vector fields, E s i m is a similarity measure between 
two images and |.| is a distance between two vector fields (L2 distance), h 1 
estimates intensity matches between T(x) and R(x). h is a trade-off between 
geometric constraints and intensity matches. (See figured]). 
The hybrid registration operates as follows : 

1. initialize h(x) and h 1 (x) to identity. 

2. minimize E\ w.r.t h 1 using the registration algorithm of 14.21 

3. minimize Ei w.r.t h using the registration algorithm of 14.31 

4. set hi = h 

5. IF CONVERGENCE REACHED STOP ELSE GOTO 2 

The intensity information is introduced in the first term of eq. (fTJ) by the 
similarity energy E s i m . 

The geometric constraint is introduced by the first term of eq. $2$ which is 
a geometric distance between model points m* = ft ( _1 (mi) (at iteration t) and 
scene points Sj . The weight Wij determines the correspondence between the two 
point sets and can be estimated by several criteria such as ICP or softAssign 
[T2] , (see figure HD 

4.2 Registration based on the minimization of E\ 

Physically-based models are often used to ensure the coherence of the deforma- 
tion field. In an elastic model, elastic forces are proportional to the displacement 
field while in a fluid model, fluid forces are proportional to the rate of change 
of the deformation field. The fluid rcgularization allows the estimation of large 
deformations while maintaining the topology of the deformation field [13j . An 
elastic model does not guarantee to keep the coherence of the deformation field. 
However fluid models are often under-constrained and the resulting deforma- 
tion field can be non-satisfactory. In our implementation we both use an elastic 
model and a fluid regularizer ( [14j ) . 

Our method implements two levels of regularization and the registration is 
done as follows : 
Do Until Convergence 
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1. Compute V p i x )Ei[h I (x) o (x +p(x))] = V c Ei. 

2. Compute a correction : c(x) = k x V C E\ * G a (x). Where k is a scalar 
controlling the magnitude of the vector field c(x) and G a (x) is a Gaussian 
kernel of standard deviation a. 

3. Compose the correction field to the current deformation. 

K( x ) = h i-i( x ) ° ( x + c ( x )) 

4. Regularize the deformation field h 1 (x) by decreasing E e i as (elastic regu- 
larization) 

where E e i as is the linear elasticity potential : 

E elas = |*±/i||V.u'|| a + f£||Vtifl| a (3) 

where h 1 (x) = x + m / (x) and u\ denotes the i th component of u (x). During 
the first step, following ([T5], [H]), we look for a small perturbation p{x) that 
minimizes E[h{x) o (x + p(x))] and call it C E\. In a second step, we compute 
a correction field by taking a fraction of the smoothed gradient V c I?i. The 
smoothing allows the filtering of noisy values of the gradient and allows to 
ensure the invertibility of the deformation field. This smoothing can be related 
to a fluid regularization controlcd by the standard deviation of the Gaussian 
kernel (a). The scalar k is tuned in order that the maximum magnitude of c(x) 
is lower than 0.5 voxel. During the step 3 the correction field c(x) is composed 
to the current deformation field. During the last step a regularization of the 
deformation field is achieved by the minimization of E e i as . This minimization 
results to a classical diffusion equation which is solved using a Gauss Seidel 
algorithm and a semi implicite scheme j!6j . The diffusion time tunes the elastic 
registration. 

4.3 Registration based on the minimization of E 2 

For the minimization of energy E2 , the algorithm uses a family of tensor product 
splines to model the deformation field. 

h(x,p) = x + ^ u km iB k (x )Bi(x 1 )B m (x 2 ) (4) 

Ukim are the spline deformation coefficients which comprise a parameter vector 
p. Bk, Bi, B m are the spline basis functions. A regularization is used to put prior 
constraints on spline parameters. The cost function to be minimized becomes: 

r 3 

C(p) = E 2 (h(x,p)) + / ]T||V U4 || 2 (5) 

J i 

where h(x) = x + u(x) and Uj(a;) denotes the i th component of u(x). The 
approximation of the integral is based on finite differences and the minimization 
of C(p) (HJ) is done by gradient descent. 
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Figure 3: Example of atlas generated using only intensity information (left) and 
using intensity and geometric information (right) 




4.4 Inversion method for computing m\ 

At each iteration, model points evolves according to the backward transforma- 
tion h~ 1 (x). We denote by m 2 ; = the initial position of surface points (at 
iteration 0). hT 1 {mi) = m\ is the transformed point at iteration t. The inversion 
is realized for each point m„ by minimizing the criterion f(jp) = \\h{p) — mi\\ 2 
with respect to p. The minima of f(p) is reached for p = h~ 1 {mi) = m*. We 
start the minimization by setting p = . 

5 Model Building 

As introduced previously, the model includes : 

• the atlas (the mean image and associated contours) 

• a shape model of the prostate embedded in the atlas 

• prior probabilities of match to user provided points (base and apex) 

5.1 Atlas construction 

Atlas-based segmentation has become a standard approach to organ delineation 
for many fields of medical image analysis particularly for the study of the brain 
(CO)- In these methods a template (atlas) is registered to a new patient image. 
The spatial transformation obtained from the registration process is used to 
map anatomical segmentations of the atlas onto the new subject. The main 
issue of these methods is the template selection. A common approach to con- 
struct an atlas is to label a particular image (Talairach atlas in neurosurgery 
for instance). Another widely used method is to construct an anatomical at- 
las from a population by averaging multiple co-registered images. Our method 
constructs an anatomical atlas by registering a population to a reference R man- 
ually chosen among the population. The reference was selected to be close to the 
mean prostate size. First, for each image Ii of the population an inhomogeneity 
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correction is performed. Then, all individuals are elastically co-registered to the 
reference R. 

After registrations, we have a set of transformations Tj that map the ref- 
erence onto each individual Ii. The atlas image is obtained by averaging the 
registered images Z; o Tj and the atlas contour (prostate) is obtained by taking 
the segmentation of the reference. Many works use an averaging of the regis- 
tered segmentations Si o T t to obtain the mean segmentation (where Si denotes 
the segmentation of the image i) or a probability map of the atlas structures [17] 
[T5] . We think that these methods are more interesting when the registration of 
the atlas is used as a prior information for another segmentation method. This 
will not improve the atlas contour in our case. The atlas construction method 
described is then iterated by taking the current atlas image as a new reference; 
we restart all registrations and the atlas construction until the mean image does 
not change significantly. 

Very often, anatomical atlas construction only makes use of intensity infor- 
mation to estimate the transformation to the reference. However, such methods 
are not usable in our case because of too important mismatches which may ap- 
pear between reference segmentation and deformed template segmentation. The 
construction of an atlas using this method gives poor results. To address this 
problem we use the hybrid method introduced in section 14.11 The geometric 
criterion (ICP) works on the contours given by the expert segmentations while 
the intensity based criterion (SSD) works on the images using a SSD as a sim- 
ilarity measure. The matching problem is then solved using an ICP criterion 
("hard assignment"): for each scene points Sj, = 1 if rrii is the closest model 
point of Sj else Wij = 0. 

5.2 Shape Model 

Our method incorporates statistical shape information in order to improve the 
quality of the resulting segmentations |19j . When computing a statistical shape 
model, a fundamental problem is the determination of a set of corresponding 
points between instances (cloud of points) of the population. In the previous 
section, we determined a set of correspondences (transformations) which maps 
the reference onto each individual, based on a framework which couples geomet- 
ric and intensity based registration. These dense transformations are now used 
to generate a set of corresponding points by deforming the cloud of points of the 
reference onto each individual. The resulting sets of corresponding points are 
statistically analyzed by principal component analysis to build a shape model. 

5.3 Region of interest (Base and Apex) 

To guide a registration algorithm, one can couple landmark and intensity-based 
registration ([H]). This is very useful in presence of large deformations to forbid 
the algorithm to fall into local minima or to incorporate a strong knowledge (ex- 
pert) when image information is poor. However the determination of landmarks 
is often very subjective. In our method instead of using landmarks we prefer to 
model the correspondence between a user-provided point and a point which be- 
longs to the surface in a probabilistic framework. We defined two interest areas 
on the prostate (base, apex) by assigning a prior probability (of correspondence 
with the user provided seed) to each atlas model point. The expert user is then 
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Figure 4: Examples of expert and automatic segmentations : top row expert 
segmentations, bottom row automatic segmentations; from the left to the right 
coronal slice, transverse slice, sagittal slice 



asked to determine one point on each region. The correspondence between the 
user provided point in the study image and the prostate surface in the atlas vol- 
ume is then carried out iteratively in a probabilistic way as described in section 

m- 

A simple method to assign prior probabilities of matches (wij ) to each couple 
(rrii,Sj) =(model point, user-provided point) is to manually define a point dj 
in a center of the region j and to assign probabilities according to the distance 
to this point using a shape function G (for example a Gaussian kernel). The 
denominator of eq ^ is used for normalization : 

_ G(rm - dj) 

~ E,: G ( m ' ~ d j) [ ' 



6 Prostate segmentation 

6.1 Anatomical Atlas To Subject Registration 

The algorithm starts with two interest points (user-provided points) and corre- 
sponding prior probability of match defined for each atlas model points. The 
first step consists in rigidly registering the atlas to the patient image using an 
intensity-based registration method using the SSD as similarity measure. In a 
second step we use the hybrid registration algorithm for the non rigid registra- 
tion of the atlas onto the study image. The similarity measure used in the hybrid 
registration is also a SSD. The model points M* = {mo, ...,tojv} represent the 
prostate surface embedded in the atlas. The scene points S = {sq, ...,sm} are 
user provided points. We recall that to* = (mi) is the moving model point 
at iteration t and itij be the prior probability of match of the scene point Sj to 
the model point m*. 
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Figure 5: Automatic segmentation - deformated atlas ; from the left to the 
right : 3D view of automatic segmentations, fusion of deformed atlas (top right 
/ bottom left) and study image (top left / bottom right) 

6.1.1 Geometric Criterion 

Here we present the method used to determine correspondence between user- 
provided (sj) points and model points (to,; ; prostate) using prior probabilities 
of match (regions of interest) defined in section (|5.3j) . The framework is based 
on the soft assign method [12| . We now consider a "binary random match 
matrix" Wij which randomly associates the scene point Sj to model point m' 
(if Wij = 1 and other elements of the row are set to zero, this means that the 
point Sj corresponds to the point m*). defined in section f5.3|) represents 
a prior probability of match P(Wij = 1) associated to each couple of points 
(m*, Sj). We now model the observation law (likelihood) as a Gaussian process, 
if a model point m* corresponds exactly to a scene point Sj , the likelihood of Sj 
is: 

P^Wa = 1) « exp (_ ll^- m iH 2 | (7) 
Weights Wij introduced in eq. ([2|) are determined by computing P(Wij = l\sj) : 

Ei ^e X pj 

Wij (9) 

These weights are recomputed at each iteration of the minimisation of Ei by 
the algorithm presented in section 14.31 

6.2 Shape Model Projection 

A shape model of the prostate has been built from a set of corresponding points 
obtained by deforming the mesh of the reference onto each individual (using 
transformations obtained during the atlas construction). At the end of the 



P(W ij = l\s j ) = 
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Table 1: Volume Based Metrics; Sensitivity (SENS) and Positive Pre- 
dictive Value (PPV) with respect to one expert segmentation 



Image 


Sens 


PPV 


Image 


Sens 


PPV 


1 


0.91 


0.75 


9 


0.91 


0.79 


2 


0.85 


0.96 


10 


0.81 


0.77 


3 


0.95 


0.6 


11 


0.94 


0.75 


4 


0.92 


0.86 


12 


0.87 


0.78 


5 


0.98 


0.73 


13 


0.89 


0.79 


6 


0.94 


0.76 


14 


0.85 


0.84 


7 


0.87 


0.47 


15 


0.87 


0.87 


8 


0.84 


0.94 


16 


0.77 


0.83 








17 


0.94 


0.72 



Table 2: Distance Based Metrics (the mean size of the prostate in 

THE AXIS BASE/ APEX IS ABOUT 48 mm) 



Apex 


Central Zone 


Base 


All 


Mean (mm) 


STD 


Mean (mm) 


STD 


Mean (mm) 


STD 


Mean (mm) 


STD 


2.91 


1.34 


2.4 


1.29 


4.30 


2.00 


3.39 


1.95 



atlas construction the reference is the atlas. It results that we have a shape 
model of the prostate surface which is embedded in the atlas. During the atlas 
to study matching the atlas mesh is deformed to segment the prostate on the 
study image. To regularize this segmentation which can be noisy, we can now 
compute the most probable shape according to the shape statistics. To this 
end, we estimate modes of variation and then we apply constraints to modes, 
to ensure plausible shapes. This is done by limiting their absolute value to be 
less than (Aj is the cigen-value of mode i) 

7 Results 

The method has been tested with the same dataset than the one used for atlas 
and shape model construction. The "leave-one-out" method has been used 
for this purpose (i.e.: the evaluated data - the patient to be automatically 
segmented - is removed from the set of data used for model construction). 
The reference image used for atlas building is not used for validation to avoid 
additionnal bias when choosing a new reference; it results that we have 17 
patients for tests. 

Results obtained using a volume-based metrics are presented in Table [U 
The two distances used are the "sensitivity" and the "positive predictive" value 
defined as follows : 

PPV- TP QTp/xr Q — TP 

TP+FP UIJJ "" TP+FN 
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Where TP denotes the number of true positives, FP denotes the number of 
false negatives and FP denotes the number of false positives. Results obtained 
using a distance based metrics are presented in Table [5] We have divided the 
prostate in 3 zones to show how the algorithm performs spatially. Results show 
that the base of the prostate is more difficult to segment. The accuracy in the 
apex zone and in the central zone is good. The convergence of the registration 
algorithm in the base of the prostate is difficult due to the large interindividual 
variability observed. For small prostates (volume less than 25 cc for a mean 
volume of 41 cc) the algorithm behaves less accurately (mean error of about 5.5 
mm in 3 cases) . We are very confident in our ability to solve this problem using 
several atlases depending on the size of the prostate in the study image. 

8 Conclusion 

This paper has presented an atlas-based semi-automatic prostate segmentation 
method for MRI. The method makes use of an hybrid registration framework 
which can deal with geometric and image-based information for both atlas con- 
struction and atlas-based segmentation. Results are very promising expect for 
small prostates (volume less than 25 cc for a mean volume of 41 cc) where the 
convergence is very difficult (mean error of about 5.5 mm in 3 cases). The use of 
an additional atlas built from a population of small prostates will certainly help 
to solve the problem. The hybrid registration allows for an efficient estimation of 
correspondences across individuals of the database. These correspondences are 
used for atlas construction and for shape model building. Many works use land- 
marks to introduce geometric constraints in intensity-based registration meth- 
ods [10] [UJ . However the frequent limitation of these approaches is the deter- 
mination of hard correspondences between reference and template image which 
is often very subjective. This work overcomes these limitations by iteratively 
estimating correspondences between the surface embedded in the atlas, which 
deforms during the registration and few user-provided points. Correspondences 
arc estimated in an elegant way using a probabilistic framework and prior infor- 
mation on expected correspondences between user-provided points and prostate 
surface. An extension of this work could be to use additional geometric features 
automatically extracted from the image by some image processing techniques 
directly in the hybrid registration framework. Current work deals with the 
addition of other expert segmentations into the database and with qualitative 
evaluation of the method from the clinical side. 

References 

[1] Jcmal, A., Sicgel, R., Ward, E., Hao, Y., Xu, J., Murray, T., Thun, M.J.: 
Cancer statistics, 2008. CA: a cancer journal for clinicians 58(2) (2008) 
71-96 

[2] C. Reynicr, J. Troccaz, P.F.A.D.C.G.J.J.D.M.B.J.Y.G.: Mri/trus data 
fusion for prostate brachytherapy. preliminary results. Medical Physics 31 
(2004) 1568-1575 



11 



[3] Vincent Daneen, J. Gastaldo, J.G.P.F.J.D.M.B.D.C.J.T.: Mri/trus data 
fusion for brachytherapy. The International Journal of Medical Robotics 
and Computer-Assisted Surgery, 2 (2006) 256-261 

[4] Zwiggclaar, R.Z.: Semi-automatic segmentation of the prostate in mri. In: 
Lcct. Notes on Comput. Sci. (2003) 

[5] Zhu, Y. Zwiggelaar, R.W.S.: Segmentations of volumetric prostate mri 
data unsing hybrid 2d + 3d shape modelling. In: Proceedings of Medical 
Image Understanding and Analysis. (2004) 

[6] Zhu. Y. Zwiggelaar, R.W.S.: An hybrid asm approach for sparse volumetric 
data segmentation. Pattern Rccogn. Image Anal 15 (2005) 346-349 

[7] Klein, S. van der Heide, U.L.I.v.V.M.S.M.P.: Automatic segmentation of 
the prostate in 3d mr images by atlas matching using localized mutual 
information. Medical Physics (2008) 

[8] Maria Jimena Costa, Herv Delingette, S.N.N. A.: Automatic segmentation 
of bladder and prostate using coupled 3d deformable models. In: Proc. 
MICCAI (2007). (2007) 

[9] HM, H.A.L.: 3d prostate boundary segmentation from ultrasound images 
using 2d active shape models. In: Conf Proc IEEE Eng Med Biol Soc. 
(2006) 

[10] P. Hellier, C.B.: Coupling dense and landmarks-based approaches for non 
rigid registration. IEEE TMI 22 (2003) 217-227 

[11] A. Azar, C. Xu, X.P.N. A.: An interactive hybrid non-rigid registration 
framework for 3d medical images. In: ISBI. (2006) 

[12] Rangarajan, A., C.H..B.: The softassign procrustes matching algorithm. 
IPMI (1997) 29-42 

[13] Thirion. J.: Image matching as a diffusion process: An analogy to maxwell 
demons. Medical Image Analysis (1998) 

[14] Radu Stefanescu, X.P., Ayache, N.: Grid powered nonlinear image regis- 
tration with locally adaptive regularization. Medical Image Analysis (2003) 

[15] Trouv: Diffcomorphisms groups and pattern matching in image analysis. 
International Journal of Computer Vision 28 (1998) 213; 221 

[16] Press, W., F.B.T.S.V.W.: Numerical Recipes in C: The Art of Scientific 
Computing. (1993) 

[17] M De Craene, A du Bois dAische, B.M.: Multi-subject registration for 
unbiased statistical atlas construction. MICCAI (2004) 

[18] H Park, PH Bland, A.H.C.M.: Least biased target selection in probabilistic 
atlas construction. LECTURE NOTES IN COMPUTER SCIENCE (2005) 

[19] Cootes, T., .H.A.: Use of active shape models for locating structures in 
medical images. Image Vision Computing (1994) 355-366 



12 



